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1.  Introduction 

There  have  been  only  a  few  parametric  models  extensively  examined 
for  application  to  reliability;  these  include  the  exponential. distri¬ 
bution  of  Epstein-Sobel  [6],  the  Weibull  distribution  [14],  and  the 
fatigue  model  of  Bimbaum-Saunders  [4],  The  one  most  widely  utilized 
1  for  electronic  components  has  been  the  exponential  model,  not  only 

because  of  its  simple  and  intuitive  properties  but  also  because  of  the 
extent  of  the  estimation  and  sampling  procedures  which  have  been  developed 

* 

from  the  theory. 

One  of  the  early  discoveries  was  that  mixtures  of  exponentially 
1  distributed  random  variables  have  a  decreasing  failure  rate,  see  [11]. 

Thus  any  two  groups  of  components  with  constant,  but  different,  failure 

rates  would,  if  mixed  and  sampled  at  random,  exhibit  a  decreasing  failure 

t 

rate.  As  a  consequence,  the  family  of  life  lengths  with  decreasing 
failure  rate  certainly  arises  in  practice  and  particular  subsets  of  this 
^  family  could  be  of  great  utility  for  specific  applications,  see  e.g. 

Cozzolino  [5].  We  examine  one  such  model  with  shape  and  scale 
parameters,  call  them  a  and  0  respectively,  which 'is  based  upon  a 
I  particular  mixture  of  exponential  distributions.  This  family  was  intro¬ 

duced  by  Af anas' ev  [1]  and  later  by  Lomax  [10]  as  a  generalization  of  a 
Pareto  distribution.  Section  3  compares  this  mixed  exponential 
distribution  to  the  exponential  distribution  using  data  from  Poseidon 
flight  control  packages. 

Kulldorff  and  VSnnman  [9]  and  Vannman  [13]  have  studied  a  variant 
of  this  mixed  exponential  model  containing  a  location  parameter.  They 
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obtained  a  best  linear  unbiased  estimate  of  the  scale  parameter 
assuming  that  the  shape  parameter,  call  it  a,  was  known  and  in  a 
region  restricted  so  that  both  the  mean  and  the  variance  exist,  namely 
a  >  2.  When  this  restriction  of  a  >  2  cannot  be  met  an  estimate  based 
on  a  few  order  statistics,  which  are  optimally  spaced,  is  claimed 
to  be  an  asymptotically  best  linear  unbiased  estimate  and  tables 
of  the  weights  as  functions  of  the  number  of  spacings  are  provided.  In 
all  cases,  the  shape  parameter  was  assumed  known  and  the  sample  was 
either  complete  or  type  II  censored.  It  is  contended  that  BLUE 
estimates  of  the  shape  parameter  are  not  attainable. 

Harris  and  Singpurwalla  [7]  examined  the  method  of  moments  as  an 
estimation  procedure  for  this  same  model  but  again  with  the  shape 
parameter  restricted  to  a  >  2  and  with  a  complete  sample. 

In  both  papers  [9]  and  [7],  it  is  stated  that  maximum  likelihood 
estimates  are  difficult  to  obtain.  In  a  later  paper  Harris  and 
Singpurwalla  [8]  exhibit  the  maximum  likelihood  equations  for  complete 
samples . 

In  this  paper  the  maximum  likelihood  estimates  for  both  the  shape 
and  scale  parameters  are  obtained,  jointly  and  separately,  with  sinple 
sufficient  conditions  given  for  their  existence.  These  estimates  are 
derived  for  censored  data  (and  a  fortori  for  complete  samples)  even  with 
a  paucity  of  failure  observations,  namely  one. 

The  existence  conditions  obtained  here  for  the  maximum  likelihood 
estimates  apply  even  to  the  case  where  the  variance  and  possibly  the 
mean  do  not  exist:  0  <  a  <2.  Moreover,  the  estimates  of  the  shape 
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parameter  a  which  have  been  obtained  from  actual  data  indicate  that 
this  region  0  <  a  <  2  is  important  because  all  the  estimates  obtained 
of  a  have  beer,  less  than  unity . 

2.  The  Model 

We  postulate  that  the  underlying  process  which  determines  the  length 
of  life  of  the  component  under  consideration  is  the  following:  The 
quality  of  construction  determines  a  level  of  resistance  to  stress  which 
the  component  can  tolerate.  The  service  environment  provides  shocks 
of  varying  magnitude  to  the  component, and  failure  takes  place  when  for 
the  first  time  the  stress  from  an  environmentally  induced  shock  exceeds 
the  strength  of  the  component. 

If  the  time  between  shocks  of  any  magnitude  is  exponentially 
distributed  with  a  mean  depending  upon  that  magnitude  then  the  life 
length  of  each  component  will  be  exponentially  distributed  with  a  failure 
rate  which  is  determined  by  the  quality  of  assent) ly.  It  follows  that 
each  component  has  a  constant  failure  rate  but  that  the  variability  in 
manufacture  and  inspection  techniques  forces  some  components  to  be 
extremely  good  while  a  few  others  are  bad  and  most  are  in  between. 

Let  be  the  life  length  of  a  component  in  such  a  service 
environment ,  with  a  constant  failure  rate  X  which  is  unknown.  The 
variability  of  manufacture  determines  various  percentages  of  the  X-values 
and  this  variability  can  be  described  by  some  distribution,  say  G. 
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Let  T  be  the  life  length  of  one  of  the  components  which  is 
selected  at  random  from  the  population  of  manufactured  components. 
We  denote  the  reliability  of  this  component  by  R  and  we  have 

R(t)  =  p[t  >  t]  for  t  >  0  . 


Let  A  be  the  random  variable  which  has  distribution  G.  We 
can  write 


R(t)  =  EaP 


LXX>  1 1  A  =  A 


■/ 


e"AtdG(A) 


(1) 


Because  of  having  a  form  which  can  fit  a  wide  variety  of  practical 
situations  when  both  scale  and  shape  parameters  are  disposable,  it  is 
assumed  that  G  is  a  garrma  distribution,  i.e.,  for  some  a  >  0,  6  >  0. 


g(M 


x<.-ie-\/ei 
r(a)  pa 


for  X  >  0 . 


That  this  assumption  is  robust,  even  when  mixing  as  few  as  five 
equally  weighted  A's,  has  been  shown  by  recent  work  of  Sunjata  in  an 
unpublished  thesis  [12].  It  follows  from  equation  (1)  that  the  reliability 
function  is 

R(t)  =  - i -  =  e"a  £n(1+tP)  •  (2) 

(l+tp)a 

The  failure  rate,  hazard  rate,  can  be  shown  to  be 

^  *  I+tp  *  (3) 

which  is  a  decreasing  function  of  t  >0. 
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Maximum  likelihood  estimates  for  a, 6  and  hence  R(t)  and  q(t)  are 
given  in  Section  5. 

3.  A  Comparison  of  the  Mixed  Exponential  with  Exponential  Using  Real  Data 

Data  has  been  accumulating  for  years  in  the  assessment  of  the 
reliability  of  electronic  equipment  for  which  there  was  no  adequate 
statistical  model.  The  following  difficulties  were  recognized  by 
practitioners:  1.  The  assumption  of  constant  or  increasing  failure  rate 
seemed  to  be  incorrect.  2.  However,  the  design  of  this  electronic 
equipment  indicated  that  individual  items  should  exhibit  a  constant 
failure  rate.  A  mixed  exponential  life  distribution  accounts  for  both 
,  the  design  knowledge  and  the  observed  life  lengths.  Maximum  likelihood 

procedures  allow  for  joint  estimation  of  the  parameters  of  this 
distribution  in  the  most  commonly  encountered  situation  where  complete 

♦  data  is  not  available. 

We  now  give  some  actual  data  sets  from  two  different  lots  of 
Poseidon  flight  control  electronic  packages  which  illustrate  these 

♦ 

points.  Each  package  has  recorded,  in  minutes,  either  a  failure  time 
or  an  alive  time.  An  alive  time  is  sometimes  called  a  "run-out"  and 
is  the  time  the  life  test  was  terminated  with  the  package  still  functioning 

First  Data  Set 

Failure  times:  1,  8,  10 

1  Alive  times:  59,  72,  76,  113,  117,  124,  145,  149,  153,  182,  320. 

Second  Data  Set 
Failure  times;  37,  53 


Alive  times: 


60,  64,  66,  70,  72,  96,  123. 
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If  we  assume  that  the  data  are  observations  from  an  exponential 
distribution  (constant  failure  rate  X)  then  using  the  total  life 
statistic j  we  have  the  estimates  of  reliability  given  in  the  left 
hand  side  of  the  table.  If  we  assume  that  the  data  are  observations 
from  the  mixed  exponential  distribution  of  equation  (2)  then  using 
estimation  techniques  derived  subsequently  in  this  paper  we  have  the 
estimates  for  reliability  given  in  the  right  hand  side  of  the  table. 

Exponential  estimate  Mixed  exponential  estimate 

of  reliability  of  reliability 


time  t 

Set  1 

Set  2 

Set  1 

Set  2 

In  min. 

\(t) 

R2(t) 

yt> 

R2(t)  • 

6 

.988 

.981 

•  915 

.976 

10 

.980 

.969 

.896 

.961 

30 

.943 

.911 

-.855 

.896 

50 

.906 

.856  . 

.836 

.  843 

100 

.821 

— 

.810 

— 

130 

.774 

— 

.801 

— 

A 

X: 

.00017 

.00312 

a:  .0453  ~ 

.  420 

p:  1.03 

.01 

Looking  at  the  data  from  the  two  sets  we  would  expect  that  at  least  for 
the  first  fifty  minutes  the  reliability  estimate  for  the  second  set  of 
data  would  be  higher  then  the  reliability  estimate  for  the  first  set  of 
data,  because  in  the  first  set  3  failure  out  of  14  trials  have  occurred 
in  the  first  ten  minutes  while  in  the  second  set  only  1  failure  out  of  9 
trials  has  occurred  in  the  first  fifty  minutes.  However,  under  the 
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exponential  assumption  the  reliability  estimates  for  the  first  data  set 
are  consistently  higher.  Mote  that  the  mixed  exponential  estimates 
are  more  consistent  with  what  the  data  show;  that  is,  for  at  least  the 
first  50  minutes  we  expect  the  reliability  estimate  for  the  second  set 
of  data  to  be  higher  than  the  reliability  estimate  for  the  first  set 
of  data.  Beyond  this  time,  however,  say  at  100  minutes,  the  data 
indicate  that  the  reliability  estimate  from  the  first  set  of  data  should 
be  higher  than  the  reliability  estimate  from  the  second  set  of  data. 

Using  mixed  exponential  estimates  this  is  the  case. 

A  statistical  test  to  determine  whether  the  data  require  a  constant 
or  decreasing  failure  rate  was  run  on  the  data  from  Sets  1  and  2.  For 
data  Set  1  we  reject  constant  failure  rate  in  favor  of  decreasing  failure 
rate  at  the  .10  level.  For  data  Set  2  we  cannot  reject  the  constant 
failure  rate  assumption.  In  this  case,  however,  the  constant  failure 
rate  estimates  for  reliability  and  the  mixed  exponential  estimates  for 
reliability  are  close.  For  data  Set  2  one  should  not  estimate 
reliability  much  beyond  about  70  minutes  since  we  do  not  have  data  to 
support  those  estimates. 

4 .  Residual  Life  Property  of  the  Model 

An  important  property  of  this  model  is  that  residual  life  on  a 
component  is  distributed  as  a  mixed  exponential.  Thus  a  ,fbum-in" 
test  of  a  component  will  yield  a  residual  life  which  is  also  in  the  same 
family.  This  property  seems  to  be  shared  only  with  the  exponential  among 
common  parametric  families  of  life  distributions. 
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The  residual  life  of  a  component  is  defined  to  be  the  life 
remaining  after  time  h,  given  that  the  component  is  alive  at  time  h. 

It  can  be  shown  that: 

A  bum- in  for  h  units  of  time  on  a  component  with  initial 

life  determined  by  a  mixed  exponential  distribution  with  parameters  . 

a  and  8  will  yield  a  residual  life  and  will  be  destributed  as  a 

mixed  exponential  with  parameters  a  and  8 

1  +Bh 

It  follows  that  this  life  length  model  is  "used'  better  than  new" 
or  "new  worse  than  used"  in  the  sense  that  we  have  stochastic 
inequality  between  a  new  component  and  one  that  has  been  burned  in, 
namely 

st 

T  <_  Th  for  all  h  >  0 . 

An  important  consequence  of  this  property  is  that  one  can  calculate 
the  value  of  the  increased  reliability  attained  by  bum- in  procedures 
as  compared  with  the  cost  of  conducting  them.  It  has  long  been  the 
practice  to  bum  in  electronic  components  based  on  intuitive  ideas  of 
"infant  mortality"  in  order  to  provide  reasonable  assurance  of  having 
detected  all  defectively  assembled  units.  This  model,  whenever  it  is 
applicable,  makes  possible  an  economic  analysis.  A  variation  of  this 
result  has  been  discussed  in  [3]. 

Example 

As  an  example  of  the  applicability  of  this  property,  consider  test 
data  from  Trident  flight  control  packages . 


Assume  that  bum-in  data  is  distributed  as  a  mixed  exponential 
with  shape  parameter  a  and  scale  parameter  8.  These  parameters  were 
estimated  ( formulas  in  Section  5  )  to  be 


a  =  .57 

8  =  .0104  2 

g  =  B 

f  1+(48)(60)B 


8  =  .0104 

After  48  hours  of  bum-in  the  residual  life  T,lQ  .  should  be  mixed 

4o  hours 

exponential  with  parameters  a  and  8  .  The  first  graph  shows 

^  '  1+28808 

the  change  in  8^.  as  a  function  of  bum-in  hours.  The  second  graph 
shows  the  ciWige  in  estimated  reliabilities  at  20  minutes  as  a  function 
of  bum-in  time,  where  reliability  at  time  20  minutes  is  estimated  to 
be 

R(t)  =  [1  +  20Bf]_  “ 

A 

Bj,  decreases  as  bum-in  time  increases. 


Data  consistent  with  success  and  failure  data  was  obtained  from  another  test 


called  Pre-Test.  We  assume  that  the  time  to  failure  T  of  flight  control 
packages  subjected  to  this  type  of  test  environment  follows  a  mixed 
exponential  distribution  with  shape  parameter  a  and  scale  parameter 
6.  Using  Pre-Test  data  these  parameters  were  estimated  (formulas  in 
Section  5  )  to  be 

a  =  .5739 
B  =  .1106 


B  =  .1106 


After  60  minutes 
exponential  with 
parameters  by 


of  test  the  residua^.  life  TgQ 

Q 

parameters  a  and  - 

1+60B 


should  be  a  mixed 
We  estimate  these 


.5739 

A 

B 

- y- 


1+tB 


.1106 

1  +  60( .1106) 


.0145 
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.0145 


a  *  .5739 
B  =  .1106 
R(l)  =  .94 

Since  R(t)  =  (1  +  Bt)-01,  we  estimate  reliability  at  1  minute  for  a  package 
which  has  not  gone  through  Pre-Test  to  be 


R£  (60)  *  [1  +  (.1106 )(!)]“• 5739  =  .94 


We  estimate  reliability  at  1  minute  for  a  package  which  has  gone  through 
Pre-Test  to  be 

Rgf(60)  =  [1  +  (.0145) (l)]"’5739  =  .99 


Note  that  these  reliabilities  are  for  Pre-Test  environments. 
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Effects  of  various  environmental  or  bum- in  tests  (in  the  sense  of 
reliability  gain)  can  be  estimated  by  comparing  reliability  estimates  forecast 
at  the  end  of  the  Job  Stack  which  proceeds  the  environment  to  reliability 
estimates  for  the  Job  Stack  which  follows  the  environment.  For  example 


/v 

a  is  fixed  at  a 


If  bum-in  is  effective  then  we  would  expect  that 


A 


< 


■■■  - 
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5.  Estimation  of  Parameters  with  Censored  Data 

Let  us  assume  throu^iout  this  section  that  we  are  given 

t, , . . . ,t  as  observed  times  of  failure  while  t,  , , , . . .  ,t  are 
i  k  k+1*  5  n 

observed  alive-times  both  obtained  from  a  mixed  exponential  (  ct,3  ) 
life  distribution  with  1  5  k  5  n.  We  define  two  functions  for  x  >  0. 


n 

Sl(X)-  =  i-i£n(1+tlX)  1  S2(X) 


1 

k 


£  (1+t.x) 
1=1  1 


A  result  on  the  maximum  likelihood  estimation  Cm.l.e.)  of  the  unknown 
parameters  is  now  given  which  utilizes  data  of  this  type. 

Theorem:  Under  the  assumptions  and  conditions  given 


(i)  When  p  >  Q  is  known,  there  exists  a  unique 

A 

m.l.e.  of  ..a  ,  say  a  ,  given  explicitly  by 

a  =  k/S.CB)  . 

(ii)  When  a  >  0  is  known,  there  exists  a  unique 

A 

m.l.e.  of  8  ,  say  8  ,  given  explicitly  by 

B  -  A1 (0) 

where  A  is  the  monotone  decreasing  function 
defined  by 

A(x)  =  kS2(x)  -  axS1(x)  for  x  >  0 


with  primes  denoting  derivatives. 


(iii)  When  a, 8  are  both  unknown,  the  m.l.e.'  of  8  , 
say  8  ,  is  given  implicitly,  when  it  exists 
positively  and  finitely,  by 

8  =  B" 1 (G) 


where  B 

B(x) 


is  the  function  defined  by 


S2(x)  s[(x) 

x  SjCx) 


for  x  >  0 


and  the  m.l.e.  of  a  ,  say  a  ,  is  given  expli 
citly  by 

a  •  k/ C  S)  . 


Theorem:  The  Inequality  for  1  <  k  <  n 

k  n  n 


k  I-*i  <  2 


1=1  j=l 


1=1 


(4) 


Is  a  sufficient  condition  which  a  (censored)  sample  from  a  mixed 
exponential  (a, 8)  distribution  must  satisfy  in  order  that  maximum 
likelihood  estimators  of  both  parameters  exist  both  positively  and  finitely 
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6.  Computational  Considerations 

The  question  which  now  arises  is:  what  kinds  of  sanples  will 
satisfy  condition  (4)?  If  k  =  n  we  see  (4)  is  equivalent  with 


from  which  we  have  the 

Remark:  A  complete  sample  of  failure  times  will  satisfy  (4)  if  the 
sample  standard  deviation  exceeds  the  sample  mean. 


It  can  be  shown  that  if  T  has  a  mixed  exponential  (ot,B)  distribution 

then 

E  [T]  =  CB(a-l)]  -1  for  a>l 

Var  [T]  *  a[B2(a-l)2(o-2)2r1  for  a>2 
Thus  the  standard  deviation  does  exceed  the  mean  for  these  values  of  the 
parameters  where  the  mean  E[T]  and  the  variance,  V[T],  exist. 


Remark:  A  sample  with  k  <  n  failure  times  and  the  remaining  n-k 
observations  truncated  at  tg  will  satisfy  (4)  if 


ba  >  ni 


i  +\~:  + 1 
V  n-k 


=  n 


2n-k 
1  n-k 


for  n  large 


where  n1  =  (t1+. . .+tk)/k  is  the  average  failure  time, 


In  the  calculation  of  B  the  equation,  C(B)  =  0,  must  be  solved  where 
C(B)  =  BS.[(B)S2(S)  or 


c  c  0  ) 


n  f  a  n  k 

y  -HjL  .  y  ln(i+t  s)  y  _a_ 

A  1+y  k  J 
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where  t, , . . .  ,t,  are  failure  times  and  t.  ......  ,t  are  censored  life 

1  k  k+1’  n 

times.  We  introduce  notation  for  the  sample  moments  as  follows: 


nr  *  E  1  V  •  «r  *  k  2  'j 

i"l  3-1 


:or  r  =  1,2,3,..., 


then  using  the  two  expansions,  valid  for  jxj  <  1  , 


.2  3 


£n(l+x)  *  x  -  +  -  ...  ,  jL-  =  1  -  x  +  x' 


and  substituting  into  C  and  simplifying  we  find,  upon  neglecting 
terms  of  third  order  in  3  ,  that 

2 

Cl+n10+n28) C?1+524 +c3 f )  ‘  [?1-C2B+?3e2^  =  0 


Multiplying  the  first  two  together  and  collecting  terms  yields 


C-j-  -  n1?1)3  -  U3  -  n2?1 


hiC?  2 
-4^0  8Z  =  o 


We  new  notice  that  the  condition  equation  (4),  can  be  written  in  the 
notation  of  (5)  as  £2  >  2n^£^. 

Thus  our  computational  procedure  to  decide  upon  the  parametric 
representation  of  the  distribution  governing  the  observations  which  have 
been  obtained  is  contained  in  the  following. 

Algorithm:  Given  t^,  .  .  .,  t^  as  failure  times  and  tk+1,.  .  .,  tR  as 
censored  times  from  a  mixed  exponential  (a,  6)  distribution 
(i)  Compute  the  sample  moments  n^,  n2,  ^3* 
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(ii)  If  c ^  <  >  assume  observations  from  a  constnat  failure  rate 

distribution  and  estimate  X  by 


X  = 


(iii)  If  >  assume  observations  are  from  a  mixed  exponential 

distribution  and  compute 


Po  * 


*”2  “  ^l**i 


^  ~  2i~i  2^  2_  ~  ^  2 


then  use  the  Kewton-Raphson  interation  procedure,  namely 
for  n  =  0  ,  1,  2,  .... 


Pn+1  =  pn  " 


n' 


C(pn) 


'>  P  =  11m  B  ,  and 


n  -*■  « 


n 


a  = 


n 


2  £n(l+t.p) 

J-l  J 


Practical  experience  indicates  that  the  Iteration  converges  very 
rapidly.  Since  the  functions  are  very  simple  a  small  programmable 
electronic  calculator,  such  as  the  HP-65,  can  be  used  to  obtain  these 
estimates.  Programs  for  the  HP-65  and  HP-97  are  available  from  the 
authors . 
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7.  Conclusion 

If  a  component  has  a  life  distribution  with  an  increasing 
failure  rate,  the  information  necessary  to  estimate  its  parameters 
must  contain  failure  times .  In  practice  this  means  that  virtually 
no  observed  failures,  within  a  fleet  of  operational  components,  provide 
little  information  with  which  to  assess  reliability. 

If  a  component  has  a  constant  failure  rate  then  both  failure 
times  and  alive  times  contribute  equally  to  its  estimation.  The 
preceeding  study  suggests  that  if  a  component  has  a  life  distribution 
with  decreasing  failure  rate  it  is  the  alive  times  within  the  data 
which  contribute  principally  to  the  estimation  of  the  parameters. 

The  problem  of  obtaining  the  usual  sampling  distributions  of 
the  maximum  likelihood  estimators  of  the  parameters  for  the  decreasing 
failure  rate  model  studied  seems  to  be  difficult  because  the  estimates 
are  only  implicitly  defined.  We  have  shown,  however,  that  when  they 
exist  the  MLE's  fora  and  8,  based  on  type  I  or  on  random  censoring, 
are  asymptotically  normally  distributed.  We  have  also  shown  that  the 
distribution  function  estimated  using  the  joint  MLE's  of  the  parameters 
is  surprisingly  closer  to  the  true  distribution  for  regions  of  interest 
in  reliability  theory,  than  is  the  estimated  distribution  function 
using  a  BLUE-k  estimate  for  the  scale  parameter  and  a  known  shape'  parameter 
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